rm(list=ls())

library(foreign)
library(tidyverse)
library(lubridate)
library(simcf)
library(margins)
library(ordinal)
library(corrplot)

# note that Zelig needs to be loaded from its archive version as it was removed from CRAN
library(Zelig)

df <- read_csv('merged_kaiser.csv')

# clean vars
df$date <- ymd(df$date)
df$state <- tolower(df$state)
df$insured <- ifelse(df$uninsured==0,1,0)
df$after_2014 <- ifelse(df$date >= as.Date('2014-01-01'),1,0)
df$strong_rep <- ifelse(df$pid5_rep == 1,1,0)
df$weak_rep <- ifelse(df$pid5_rep == 0.75,1,0)
df$weak_dem <- ifelse(df$pid5_rep == 0.25,1,0)
df$strong_dem <- ifelse(df$pid5_rep == 0,1,0)
df$excellent_health <- ifelse(df$health_poor == 1,1,0)
df$verygood_health <- ifelse(df$health_poor == 0.75,1,0)
df$good_health <- ifelse(df$health_poor == 0.5,1,0)
df$fair_health <- ifelse(df$health_poor == 0.25,1,0)
df$race_other <- ifelse(df$white == 0 & df$black == 0, 1, 0)

# Correlation analysis
lm.out <- lm(obamacare_fav ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), df[df$over_65 == 0,])
lm.out2 <- lm(obamacare_fav ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), df[df$over_65 == 0 & df$date < as.Date('2014-01-01'),])
lm.out3 <- lm(obamacare_fav ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), df[df$over_65 == 0 & df$date >= as.Date('2014-01-01'),])

# Appendix Table SI3.1
stargazer(lm.out, lm.out2, lm.out3,omit=c('state','survey'),
          type='html', style='APSR',
          column.labels =  c('Full','Pre-1/1/2014','Post-1/1/2014'),
          dep.var.labels=c('Oppose ACA'),
          add.lines = list(c('State FE', 'Y','Y','Y'),
                           c('Date FE', 'Y','Y','Y')),
          covariate.labels = c('Poor Health','PID 5 (R)','Ideo (Conservative)',
                               'Income Over 50 (ref= <50)','Income Missing','Age',
                               'College','Female','Black (ref=White)','Race Other','Married'),
          star.cutoffs = c(0.05, 0.01, 0.001),
          out='r1.html')

# Robustness using CLM for ordered probit Table SI3.2
df$ocare_fav_fac <- as.factor(df$obamacare_fav)
lm.out <- polr(ocare_fav_fac ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), data=df[df$over_65 == 0,], Hess=T)
lm.out2 <- polr(ocare_fav_fac ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), data=df[df$over_65 == 0 & df$date < as.Date('2014-01-01'),], Hess=T)
lm.out3 <- polr(ocare_fav_fac ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), data=df[df$over_65 == 0 & df$date >= as.Date('2014-01-01'),], Hess=T)

stargazer(lm.out, lm.out2, lm.out3,
          omit=c('state','survey'),
          type='html', style='APSR',
          column.labels =  c('Full','Pre-01/01/2014','Post-01/01/2014'),
          dep.var.labels=c('Oppose ACA'),
          add.lines = list(c('State FE', 'Y','Y','Y'),
                           c('Date FE', 'Y', 'Y', 'Y')),
          covariate.labels = c('Poor Health','PID 5 (R)','Ideo (Conservative)',
                               'Income Over 50 (ref=<50)','Income Missing','Age',
                               'College','Female','Black','Race Other','Married'),
          star.cutoffs = c(0.05, 0.01, 0.001),
          out='r1_robust.html')

# Appendix Table SI3.3
lm.outb <- lm(obamacare_fav ~ excellent_health + verygood_health + good_health + fair_health + 
                strong_rep + weak_rep + weak_dem + strong_dem +  + ideo_cons + income_over50 + 
                income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), df[df$over_65 == 0,])
lm.outb2 <- lm(obamacare_fav ~ excellent_health + verygood_health + good_health + fair_health + 
                 strong_rep + weak_rep + weak_dem + strong_dem +  + ideo_cons + income_over50 + 
                 income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), df[df$over_65 == 0 & df$date < as.Date('2014-01-01'),])
lm.outb3 <- lm(obamacare_fav ~ excellent_health + verygood_health + good_health + fair_health + 
                 strong_rep + weak_rep + weak_dem + strong_dem +  + ideo_cons + income_over50 + 
                 income_missing + age + college + female + black + race_other + married + as.factor(state) + as.factor(survey), df[df$over_65 == 0 & df$date >= as.Date('2014-01-01'),])

stargazer(lm.outb, lm.outb2, lm.outb3,omit=c('state','survey'),
          type='html', style='APSR',
          column.labels =  c('Full','Pre-01/01/2014','Post-01/01/2014'),
          dep.var.labels=c('Oppose ACA'),
          add.lines = list(c('State FE', 'Y','Y','Y'),
                           c('Date FE', 'Y','Y','Y')),
          covariate.labels = c('Excellent Health (ref=poor)','Very Good Health','Good Health','Fair Health',
                               'Strong R (ref=Ind)','Weak R','Weak Dem','Strong Dem','Ideo (Conservative)',
                               'Income Over 50 (ref= <50)','Income Missing','Age',
                               'College','Female','Black (ref=White)','Race Other','Married'),
          star.cutoffs = c(0.05, 0.01, 0.001),
          out='r1b.html')

df$survey <- as.factor(df$survey)
df$state <- as.factor(df$state)
form <- obamacare_fav ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + survey + state

runZel <- function(data, model){
  out <- zelig(form, data = data, model = 'ls') %>%
    setx(health_poor = 0) %>%
    setx1(health_poor = 1) %>%
    sim() 
  health <- data.frame(x='Poor Health',y=mean(out$sim.out$x1$fd[[1]]),
             lower=quantile(out$sim.out$x1$fd[[1]],0.025),
             upper=quantile(out$sim.out$x1$fd[[1]],0.975),
             Model=model)
  out <- zelig(form, data = data, model = 'ls') %>%
    setx(pid5_rep = 0) %>%
    setx1(pid5_rep = 1) %>%
    sim() 
  party <- data.frame(x='Party Rep',y=mean(out$sim.out$x1$fd[[1]]),
                       lower=quantile(out$sim.out$x1$fd[[1]],0.025),
                       upper=quantile(out$sim.out$x1$fd[[1]],0.975),
                      Model=model)
  
  out <- zelig(form, data = data, model = 'ls') %>%
    setx(ideo_cons = 0) %>%
    setx1(ideo_cons = 1) %>%
    sim() 
  conservative <- data.frame(x='Ideo Cons',y=mean(out$sim.out$x1$fd[[1]]),
                      lower=quantile(out$sim.out$x1$fd[[1]],0.025),
                      upper=quantile(out$sim.out$x1$fd[[1]],0.975),
                      Model=model)
  res <- bind_rows(health, party, conservative)
  return(res)
}

all <- runZel(df[df$over_65 == 0,], model='Full')
before <- runZel(df[df$over_65 == 0 & df$date < as.Date('2014-01-01'),], model='Before Rollout')
after <- runZel(df[df$over_65 == 0 & df$date >= as.Date('2014-01-01'),], model='After Rollout')

datafp <- bind_rows(all, before, after) %>%
  mutate(Model=factor(Model, levels=c('Full','Before Rollout','After Rollout')),
         x=case_when(
           x=='Ideo Cons' ~ 'Ideology \n(Conservative)',
           x=='Party Rep' ~ 'Party \n(Republican)',
           x=='Poor Health' ~ 'Poor Health'
         ),
         x=factor(x,levels=c('Poor Health','Party \n(Republican)','Ideology \n(Conservative)'))) 

datafp %>%
  ggplot(aes(x, y, shape=Model)) +
  geom_errorbar(aes(x, ymin=lower, ymax=upper), width=0.05,position=position_dodge(width=0.2)) +
  geom_point(position=position_dodge(width=0.2), fill='white',size=2.5) +
  #coord_flip() +
  theme_bw() +
  geom_hline(yintercept = 0, linetype=2, color='grey') +
  labs(x='', y='First Difference (min to max)') +
  scale_shape_manual(values=c(21,22,24)) +
  ggsave(width=5,height=4,file='plot1d.png')

# Alternate specification Self Interest

df2 <- read_csv('kaiser_2016_robust.csv')
df2$obamacare_disapprove <- (df2$obamacare_disapprove - 1)/3
df2$race_other <- ifelse(df2$white == 0 & df2$black == 0,1,0)

form <- obamacare_disapprove ~ 
  household_preexisting  + 
  medical_difficulty_cost_scale + 
  medicaid + private_insurance + 
  pid5_republican + ideology_conservative + 
 income_over50 + income_missing + 
  black + race_other + age + female + college + married + as.factor(state)
form2 <- obamacare_disapprove ~ 
  I(income <= 3)*(household_preexisting  + 
  medical_difficulty_cost_scale + 
  medicaid + private_insurance) + 
  pid5_republican + ideology_conservative + 
  black + race_other + age + female + college + married+ as.factor(state)

lm.out1 <- lm(form, data=df2[df2$over_65 == 0,])
lm.out2 <- lm(form2, data=df2[df2$over_65 == 0,])

stargazer(lm.out1, lm.out2,
          omit=c('state'),
          type='html', style='APSR',
          #column.labels =  c('Full','Lower Tercile Income'),
          dep.var.labels=c('Oppose ACA'),
          add.lines = list(c('State FE', 'Y','Y')),
          #covariate.labels = c('Pre-existing Conditions','Medical Cost Scale','Medical Anxiety Scale',
          #                     'On Medicaid','Private Insurance',
          #                     'PID 5 (R)','Ideo (Conservative)',
          #                     'Income Over 50','Income Missing','White',
          #                     'Black','Age','Female','College','Married'),
          star.cutoffs = c(0.05, 0.01, 0.001),
          out='r2_altnerative_spec_selfinterest.html')


# Cor Plot
M <- cor(cbind(df2$obamacare_disapprove, df2$health_poor, df2$medical_difficulty_cost_scale, df2$household_preexisting, df2$medicaid, df2$private_insurance, df2$pid5_republican, df2$ideology_conservative, df2$income_over50, df2$income_missing, df2$age, df2$college, df2$female, df2$white, df2$black, df2$married), use='complete.obs',)
rownames(M) <- c('Oppose ACA','Poor H','Pre-Existing Cond', 'Medical Cost Difficulty','Medicaid','Private Insurance', 'PID','Ideo Cons','Income 50','Income Miss','Age','College','Female','White','Black','Married')
colnames(M) <- c('Oppose ACA','Poor H','Pre-Existing Cond', 'Medical Cost Difficulty','Medicaid','Private Insurance', 'PID','Ideo Cons','Income 50','Income Miss','Age','College','Female','White','Black','Married')
postscript("corrplot.eps", height=12, width=12)
corrplot(M, type = 'upper',diag = F,addCoef.col = "black",tl.srt=45)
dev.off()

# Exchange Enrollment
# form <- exchange_eligible ~ health_poor + pid5_rep + ideo_cons + income_over50 + income_missing + age + college + female + black + race_other + married + state + survey
# runZel <- function(data, model){
#   out <- zelig(form, data = data, model = 'ls') %>%
#     setx(health_poor = 0) %>%
#     setx1(health_poor = 1) %>%
#     sim() 
#   health <- data.frame(x='Poor Health',y=mean(out$sim.out$x1$fd[[1]]),
#                        lower=quantile(out$sim.out$x1$fd[[1]],0.025),
#                        upper=quantile(out$sim.out$x1$fd[[1]],0.975),
#                        Model=model)
#   out <- zelig(form, data = data, model = 'ls') %>%
#     setx(pid5_rep = 0) %>%
#     setx1(pid5_rep = 1) %>%
#     sim() 
#   party <- data.frame(x='Party Rep',y=mean(out$sim.out$x1$fd[[1]]),
#                       lower=quantile(out$sim.out$x1$fd[[1]],0.025),
#                       upper=quantile(out$sim.out$x1$fd[[1]],0.975),
#                       Model=model)
#   
#   out <- zelig(form, data = data, model = 'ls') %>%
#     setx(ideo_cons = 0) %>%
#     setx1(ideo_cons = 1) %>%
#     sim() 
#   conservative <- data.frame(x='Ideo Cons',y=mean(out$sim.out$x1$fd[[1]]),
#                              lower=quantile(out$sim.out$x1$fd[[1]],0.025),
#                              upper=quantile(out$sim.out$x1$fd[[1]],0.975),
#                              Model=model)
#   res <- bind_rows(health, party, conservative)
#   return(res)
# }
# all <- runZel(df[df$over_65 == 0,], model='Full')
# 
# g.out <- all %>%
#   mutate(x=case_when(
#            x=='Ideo Cons' ~ 'Ideology \n(Conservative)',
#            x=='Party Rep' ~ 'Party \n(Republican)',
#            x=='Poor Health' ~ 'Poor Health'
#          ),
#          x=factor(x,levels=c('Poor Health','Party \n(Republican)','Ideology \n(Conservative)'))) %>%
#   ggplot(aes(x, y, shape=Model)) +
#   geom_errorbar(aes(x, ymin=lower, ymax=upper), width=0.05,position=position_dodge(width=0.2)) +
#   geom_point(position=position_dodge(width=0.2), fill='white',size=2.5) +
#   #coord_flip() +
#   theme_bw() +
#   geom_hline(yintercept = 0, linetype=2, color='grey') +
#   labs(x='', y='First Difference (min-max)') +
#   scale_shape_manual(values=c(21,22,24))
# g.out


# ggsave(g.out, width=5, height=4, file='figures1.eps')